Imports System.Drawing
<CLSCompliant(False)> Public Class ColorUtil
    'Author: Yves Vander Haeghen (Yves.VanderHaeghen@UGent.be)
    'Version: 1.0
    'VersionDate": 13 june 2003

    'MS surely fucked up its implementation of color and color transforms!!!!
    'Some methods are not static and need the class to be instantiated. 

    'Enums used as array indices to make code more readable

    'January 2004. Clarification for LUTs.
    'Integer LUTS always work from integer to integer, and remain in the same color space.
    'Float LUTs always work from integer to single, and always to linear light RGB. This means
    'they incorporate the inverse gamma correction if necessary.

    Public Enum RGB
        red = 0
        green = 1
        blue = 2
    End Enum

    Public Enum CIEXYZ
        X = 0
        Y = 1
        Z = 2
    End Enum

    Public Enum CIELab
        L = 0
        a = 1
        b = 2
    End Enum

    Public Enum WhitePoint
        CIED65 = 0
        CIED55 = 1
        CIEA = 2
        CIEC = 3
        CIED93 = 4
        CIED50 = 5
        HDTVD65 = 6
    End Enum

    Public Enum ColorSpace
        GenericRGB = 0
        GenericGammaRGB
        sRGB
        sGammaRGB
        CIEXYZ
        CIELab
    End Enum

    Public Enum ColorSpaceTransform
        'They use the number of terms as ID, so we can decude them from the 
        'transformation matrix by counting the rows
        Linear = 3
        NonLinear6Term = 6
        NonLinear8Term = 8
        NonLinear9Term = 9
        NonLinear11Term = 11
        NonLinear14Term = 14
        NonLinear20Term = 20
    End Enum

    Public Enum FloatLUTType
        RGBToRGB = 0
        GammaRGBToRGB
    End Enum

    Public Enum IntegerLUTType
        RGBToRGB = 0
        GammaRGBToGammaRGB
    End Enum

    'Rec709 transforms and white point matrices
    'Class variables
    Private sWhite(6, 2), sXYZToSRGB(2, 2), sSRGBtoXYZ(2, 2) As Single

    'Some shared variables
    Private Shared sColorSpaceTransformDesignMatrix(,), sColorSpaceTransformAimCIELabColors(,) As Single
    Private Shared iColorSpaceTransformMetricPower As Integer = 1

    Public Sub New()
        'Object constructor
        MyBase.new()

        'Initalise the LUTs and whitepoint arrays

        'Whitepoints.
        sWhite(WhitePoint.HDTVD65, CIEXYZ.X) = 94.825
        sWhite(WhitePoint.HDTVD65, CIEXYZ.Y) = 100.0# 'Poynton
        sWhite(WhitePoint.HDTVD65, CIEXYZ.Z) = 107.391

        sWhite(WhitePoint.CIED65, CIEXYZ.X) = 95.017
        sWhite(WhitePoint.CIED65, CIEXYZ.Y) = 100.0# 'Wyszecki
        sWhite(WhitePoint.CIED65, CIEXYZ.Z) = 108.813

        sWhite(WhitePoint.CIED55, CIEXYZ.X) = 95.642
        sWhite(WhitePoint.CIED55, CIEXYZ.Y) = 100.0# 'Wyszecki
        sWhite(WhitePoint.CIED55, CIEXYZ.Z) = 92.085

        sWhite(WhitePoint.CIEA, CIEXYZ.X) = 109.828
        sWhite(WhitePoint.CIEA, CIEXYZ.Y) = 100.0# 'Wyszecki
        sWhite(WhitePoint.CIEA, CIEXYZ.Z) = 35.547

        sWhite(WhitePoint.CIEC, CIEXYZ.X) = 98.041
        sWhite(WhitePoint.CIEC, CIEXYZ.Y) = 100.0# 'Wyszecki
        sWhite(WhitePoint.CIEC, CIEXYZ.Z) = 118.103

        sWhite(WhitePoint.CIED50, CIEXYZ.X) = 96.3963
        sWhite(WhitePoint.CIED50, CIEXYZ.Y) = 100.0# 'Wyszecki
        sWhite(WhitePoint.CIED50, CIEXYZ.Z) = 82.4145

        'Transform 3x3 from XYZ to sRGB using Rec 709 primaries
        sXYZToSRGB(0, 0) = 3.240479 / 100.0#
        sXYZToSRGB(1, 0) = -0.969256 / 100.0#
        sXYZToSRGB(2, 0) = 0.055648 / 100.0#
        sXYZToSRGB(0, 1) = -1.53715 / 100.0#
        sXYZToSRGB(1, 1) = 1.875992 / 100.0#
        sXYZToSRGB(2, 1) = -0.204043 / 100.0#
        sXYZToSRGB(0, 2) = -0.498535 / 100.0#
        sXYZToSRGB(1, 2) = 0.041556 / 100.0#
        sXYZToSRGB(2, 2) = 1.057311 / 100.0#

        'Transform 3x3 from sRGB to XYZ using Rec 709 primaries
        sSRGBtoXYZ(0, 0) = 0.412453 * 100.0#
        sSRGBtoXYZ(1, 0) = 0.212671 * 100.0#
        sSRGBtoXYZ(2, 0) = 0.019334 * 100.0#
        sSRGBtoXYZ(0, 1) = 0.35758 * 100.0#
        sSRGBtoXYZ(1, 1) = 0.71516 * 100.0#
        sSRGBtoXYZ(2, 1) = 0.119193 * 100.0#
        sSRGBtoXYZ(0, 2) = 0.180423 * 100.0#
        sSRGBtoXYZ(1, 2) = 0.072169 * 100.0#
        sSRGBtoXYZ(2, 2) = 0.950227 * 100.0#
    End Sub

    Public Shared Function RGBName(ByVal i As Integer) As String
        'Return the name associated with a color component
        Select Case i
            Case RGB.blue
                Return "Blue"
            Case RGB.green
                Return "Green"
            Case RGB.red
                Return "Red"
        End Select
        Return Nothing
    End Function

    Public Shared Function CIEXYZName(ByVal i As Integer) As String
        'Return the name associated with a color  component 
        Select Case i
            Case CIEXYZ.X
                Return "X"
            Case CIEXYZ.Y
                Return "Y"
            Case CIEXYZ.Z
                Return "Z"
        End Select
        Return Nothing
    End Function

    Public Overloads Shared Sub RGBTosRGB(ByVal sRGB() As Single, _
                                          ByVal sSRGB() As Single, _
                                          ByVal sTransform(,) As Single)
        'Transform linear RGB using a polynomial mapping technique.
        'RGB both in [0.0 1.0] range. Convenience function, it just makes code more readable.
        TransformTristimulus(sRGB, sSRGB, sTransform)
    End Sub

    Public Overloads Shared Sub RGBTosRGB(ByVal sRGB(,) As Single, _
                                          ByVal sSRGB(,) As Single, _
                                          ByVal sTransform(,) As Single)
        'Transform linear RGB using a polynomial mapping technique.
        'RGB both in [0.0 1.0] range. Convenience function, it just makes code more readable.
        TransformTristimulus(sRGB, sSRGB, sTransform)
    End Sub

    Public Overloads Shared Sub TransformTristimulus(ByVal sC1(,) As Single, _
                                               ByVal sC2(,) As Single, _
                                               ByVal sTransform(,) As Single)
        Dim sC1Row(2), sC2Row(2) As Single
        Dim iHiRow, i As Integer
        iHiRow = sC1.GetUpperBound(0)
        For i = 0 To iHiRow
            Algebra.GetMatrixRow(sC1, sC1Row, i)
            TransformTristimulus(sC1Row, sC2Row, sTransform)
            Algebra.SetMatrixRow(sC2, sC2Row, i)
        Next
    End Sub

    Public Overloads Shared Sub TransformTristimulus(ByVal sC1() As Single, _
                                               ByVal sC2() As Single, _
                                               ByVal sTransform(,) As Single)
        'Transform linear tristimulus sC1 to sC2 using a polynomial mapping technique.
        'Both tristimuli in [0.0 1.0] range
        ' Comment Oct. 2003: Profiling shows a lot of waste in ColorTripletToNTuple, which is
        ' not needed when dealing with a linear transform! Testing for this avoids an 
        ' unneeded array copy
        Dim iNrTerms As Integer = sTransform.GetUpperBound(1) + 1
        Dim sCNL(iNrTerms - 1) As Single

        If iNrTerms > ColorSpaceTransform.Linear Then
            'Compute color N-tuple with may contain higher-order terms.
            ColorTripletToNTuple(sC1, sCNL, iNrTerms)

            'Left multiply color N-tuple with transform
            Algebra.Product(sTransform, sCNL, sC2)
        Else
            'Left multiply color N-tuple with transform
            Algebra.Product(sTransform, sC1, sC2)
        End If
    End Sub

    Public Overloads Sub XYZToCIELab(ByVal sXYZ() As Single, _
                                     ByVal sLab() As Single)
        XYZToCIELab(sXYZ, sLab, WhitePoint.CIED65)
    End Sub

    Public Overloads Sub XYZToCIELab(ByVal sXYZ() As Single, _
                           ByVal sLab() As Single, _
                           ByVal iWhite As Integer)
        'XYZ to CIELab transform.
        Dim sY As Single

        sY = sXYZ(CIEXYZ.Y) / sWhite(iWhite, CIEXYZ.Y)
        sLab(CIELab.L) = CSng(116 * CIEf(sY) - 16.0#)
        sLab(CIELab.a) = 500 * (CIEf(sXYZ(CIEXYZ.X) / sWhite(iWhite, CIEXYZ.X)) - CIEf(sY))
        sLab(CIELab.b) = 200 * (CIEf(sY) - CIEf(sXYZ(CIEXYZ.Z) / sWhite(iWhite, CIEXYZ.Z)))
    End Sub

    Public Overloads Sub CIELabToXYZ(ByVal sLab() As Single, _
                                      ByVal sXYZ() As Single)
        CIELabToXYZ(sLab, sXYZ, WhitePoint.CIED65)
    End Sub

    Public Overloads Sub CIELabToXYZ(ByVal sLab() As Single, _
                           ByVal sXYZ() As Single, _
                           ByVal iWhite As Integer)
        'CIELab to XYZ transform.
        Dim sTemp As Single, i As Integer

        sTemp = CSng((sLab(0) + 16.0#) / 116.0#)
        sXYZ(CIEXYZ.Y) = CIEfInverse(sTemp) * sWhite(iWhite, CIEXYZ.Y)
        sXYZ(CIEXYZ.X) = CIEfInverse(CSng(sLab(1) / 500.0# + sTemp)) * sWhite(iWhite, CIEXYZ.X)
        sXYZ(CIEXYZ.Z) = CIEfInverse(CSng(-sLab(2) / 200.0# + sTemp)) * sWhite(iWhite, CIEXYZ.Z)
        For i = 0 To 2
            If sXYZ(i) < 0 Then sXYZ(i) = 0
        Next i
    End Sub

    Public Shared Function CIELabTodE(ByVal sLab1() As Single, _
                                      ByVal sLab2() As Single) As Single
        Dim sDeltaLab(2) As Single
        Algebra.Subtract(sLab1, sLab2, sDeltaLab)
        Return Algebra.Norm(sDeltaLab)
    End Function

    Public Shared Function CIELabTodC(ByVal sLab1() As Single, _
                               ByVal sLab2() As Single) As Single
        'Use only chromatic info, not luminance
        CIELabTodC = CSng(((sLab1(CIELab.a) - sLab2(CIELab.a)) ^ 2 + (sLab1(CIELab.b) - sLab2(CIELab.b)) ^ 2) ^ 0.5)
    End Function

    Public Sub CIELabTosRGB(ByVal sLab() As Single, _
                            ByVal sRGB() As Single)
        'Input must have D65 whitepoint!
        Dim sXYZ(2) As Single
        CIELabToXYZ(sLab, sXYZ, WhitePoint.CIED65)
        XYZToSRGB(sXYZ, sRGB)
    End Sub

    Public Sub sRGBToCIELab(ByVal sRGB() As Single, _
                            ByVal sLab() As Single)
        'Output has D65 whitepoint!
        Dim sXYZ(2) As Single
        sRGBtoXYZ(sRGB, sXYZ)
        XYZToCIELab(sXYZ, sLab, WhitePoint.CIED65)
    End Sub

    Public Overloads Function sRGBTodE(ByVal sSRGB1() As Single, ByVal sSRGB2() As Single) As Single
        'Output has D65 whitepoint!
        Dim sLab1(2), sLab2(2) As Single
        sRGBToCIELab(sSRGB1, sLab1)
        sRGBToCIELab(sSRGB2, sLab2)
        Return CIELabTodE(sLab1, sLab2)
    End Function

    Public Overloads Sub sRGBTodE(ByVal sSRGB1(,) As Single, ByVal sSRGB2(,) As Single, ByVal sdE() As Single)
        'Output has D65 whitepoint!
        Dim i As Integer, iMax As Integer = sSRGB1.GetUpperBound(0)

        For i = 0 To iMax
            Dim sRGB1(2), sRGB2(2) As Single

            Algebra.GetMatrixRow(sSRGB1, sRGB1, i)
            Algebra.GetMatrixRow(sSRGB2, sRGB2, i)
            sdE(i) = sRGBTodE(sRGB1, sRGB2)
        Next
    End Sub

    Public Overloads Function XYZTodE(ByVal sXYZ1() As Single, _
                          ByVal sXYZ2() As Single) As Single
        Return XYZTodE(sXYZ1, sXYZ2, WhitePoint.CIED65)
    End Function

    Public Overloads Function XYZTodE(ByVal sXYZ1() As Single, _
                            ByVal sXYZ2() As Single, _
                            ByVal iWhite As Integer) As Single
        'Compute dE from two XYZ tristimulus values. We take some shortcuts here for speed
        Dim sDLab(2) As Single

        sDLab(CIELab.L) = 116 * (CIEf(sXYZ1(CIEXYZ.Y) / sWhite(iWhite, CIEXYZ.Y)) - _
          CIEf(sXYZ2(CIEXYZ.Y) / sWhite(iWhite, CIEXYZ.Y)))
        sDLab(CIELab.a) = 500 * (CIEf(sXYZ1(CIEXYZ.X) / sWhite(iWhite, CIEXYZ.X)) - _
          CIEf(sXYZ1(CIEXYZ.Y) / sWhite(iWhite, CIEXYZ.Y)) - _
          CIEf(sXYZ2(CIEXYZ.X) / sWhite(iWhite, CIEXYZ.X)) + _
          CIEf(sXYZ2(CIEXYZ.Y) / sWhite(iWhite, CIEXYZ.Y)))
        sDLab(CIELab.b) = 200 * (CIEf(sXYZ1(CIEXYZ.Y) / sWhite(iWhite, CIEXYZ.Y)) - _
          CIEf(sXYZ1(CIEXYZ.Z) / sWhite(iWhite, CIEXYZ.Z)) - _
          CIEf(sXYZ2(CIEXYZ.Y) / sWhite(iWhite, CIEXYZ.Y)) + _
          CIEf(sXYZ2(CIEXYZ.Z) / sWhite(iWhite, CIEXYZ.Z)))
        Return Algebra.Norm(sDLab)
    End Function

    Private Function CIEf(ByVal sValue As Single) As Single
        'CIE f function for the computation of CIE Lab values

        If sValue >= 0.008856 Then
            CIEf = CSng(sValue ^ 0.3333333333)
        Else
            CIEf = CSng(7.787 * sValue + 0.1379310344828)
        End If
    End Function

    Private Function CIEfInverse(ByVal sValue As Single) As Single
        'CIE f function for the computation of CIE Lab values

        If sValue >= 0.206892706482 Then
            CIEfInverse = CSng(sValue ^ 3)
            'Debug.Print "CIE f inverse of " & sValue & ": " & CIEfInverse
        Else
            CIEfInverse = CSng((sValue - 0.1379310344828) / 7.787)
            'Debug.Print "CIE f inverse (linear part) of " & sValue & ": " & CIEfInverse
        End If
    End Function

    Public Sub XYZToSRGB(ByVal sXYZ() As Single, ByVal sRGB() As Single)
        'Transform CIE XYZ to RGB using the Rec 709 primaries (sRGB space)
        'RGB in [0.0 1.0] range, XYZ in [0.0 1xx.x] range
        Algebra.Product(sXYZToSRGB, sXYZ, sRGB)
    End Sub

    Public Sub sRGBtoXYZ(ByVal sRGB() As Single, ByVal sXYZ() As Single)
        'Transform RGB to XYZ using the Rec 709 primaries (sRGB space)
        'RGB in [0.0 1.0] range, XYZ in [0.0 1xx.X] range
        Algebra.Product(sSRGBtoXYZ, sRGB, sXYZ)
    End Sub

    Public Sub XYZToxyY(ByVal sXYZ() As Single, ByVal sXYY() As Single)
        Dim sSum As Single

        sSum = sXYZ(CIEXYZ.X) + sXYZ(CIEXYZ.Y) + sXYZ(CIEXYZ.Z)
        sXYY(0) = sXYZ(CIEXYZ.X) / sSum
        sXYY(1) = sXYZ(CIEXYZ.Y) / sSum
        sXYY(2) = sXYZ(CIEXYZ.Y)
    End Sub

    Public Shared Function GammaCorrection(ByVal sIn As Single) As Single
        'Apply gamma correction from input range [0.0 1.0] to [0.0 1.0]
        'Do clipping
        'If sIn >= 1.0# Then
        '  GammaCorrection = 1.0#
        'ElseIf sIn < 0.0# Then
        '  GammaCorrection = 0.0#
        'ElseIf sIn <= 0.018 Then
        '  GammaCorrection = 4.5 * sIn
        'Else
        '  GammaCorrection = 1.099 * sIn ^ 0.45 - 0.099
        'End If

        ' Still do clipping, but assume most will be in range and use that to speed things up
        ' This version has 30% speed gain!
        If sIn > 0.018 Then
            If sIn >= 1.0# Then
                Return 1.0#
            Else
                Return CSng(1.099 * (sIn ^ 0.45) - 0.099)
            End If
        Else
            If sIn < 0.0# Then
                Return 0.0#
            Else
                Return CSng(4.5 * sIn)
            End If
        End If
    End Function

    Public Shared Function InverseGammaCorrection(ByVal sIn As Single) As Single
        'Apply gamma correction from input range [0.0 1.0] to [0.0 1.0]
        If sIn <= 0.081 Then
            Return CSng(sIn / 4.5)
        Else
            Return CSng(((sIn + 0.099) / 1.099) ^ 2.22222222)
        End If
    End Function

    Public Shared Sub Windowing(ByVal sRGB() As Single, _
                                ByVal sOffsetRGB() As Single, _
                                ByVal sScaleRGB() As Single, _
                                ByVal sWindowedRGB() As Single)
        'Apply subwindow selection to a RGB value: sWRGB = (sRGB - sOffsetRGB) .* sScaleRGB
        sWindowedRGB(RGB.red) = Algebra.Scale(sRGB(RGB.red), sOffsetRGB(RGB.red), sScaleRGB(RGB.red))
        sWindowedRGB(RGB.green) = Algebra.Scale(sRGB(RGB.green), sOffsetRGB(RGB.green), sScaleRGB(RGB.green))
        sWindowedRGB(RGB.blue) = Algebra.Scale(sRGB(RGB.blue), sOffsetRGB(RGB.blue), sScaleRGB(RGB.blue))
    End Sub

    Public Shared Sub Windowing2(ByVal sRGBIn() As Single, _
                              ByVal sLowRGBIn() As Single, _
                              ByVal sHighRGBIn() As Single, _
                              ByVal sRGBOut() As Single)
        'Scale from given range to 0.0 - 1.0
        'Use only for small groups of values, otherwise precompute the scale and use Windowing
        Dim sLowRGBOut() As Single = {0, 0, 0}
        Dim sHighRGBOut() As Single = {1, 1, 1}
        Windowing2(sRGBIn, sLowRGBIn, sHighRGBIn, sLowRGBOut, sHighRGBOut, sRGBOut)
    End Sub

    Public Shared Sub Windowing2(ByVal sRGBIn() As Single, _
                            ByVal sLowRGBIn() As Single, _
                            ByVal sHighRGBIn() As Single, _
                            ByVal sLowRGBOut() As Single, _
                            ByVal sHighRGBOut() As Single, _
                            ByVal sRGBOut() As Single)
        'Apply scaling from one range to another.
        Dim i As Integer, temp As Single
        For i = ColorUtil.RGB.red To ColorUtil.RGB.blue
            temp = (sHighRGBOut(i) - sLowRGBOut(i)) / (sHighRGBIn(i) - sLowRGBIn(i))
            sRGBOut(i) = temp * (sRGBIn(i) - sLowRGBIn(i)) + sLowRGBOut(i)
        Next
    End Sub

    Public Overloads Shared Sub RGBToGammaRGB(ByVal sRGB() As Single, ByVal iGammaRGB() As Integer)
        RGBToGammaRGB(sRGB, iGammaRGB, 255)
    End Sub

    Public Overloads Shared Sub RGBToGammaRGB(ByVal sRGB() As Single, ByVal iGammaRGB() As Integer, ByVal iScale As Integer)
        'Gamma correction. Input RGB in [0.0 1.0], output RGB [0 iScale].
        iGammaRGB(RGB.red) = CInt(Math.Round(GammaCorrection(sRGB(RGB.red)) * iScale))
        iGammaRGB(RGB.green) = CInt(Math.Round(GammaCorrection(sRGB(RGB.green)) * iScale))
        iGammaRGB(RGB.blue) = CInt(Math.Round(GammaCorrection(sRGB(RGB.blue)) * iScale))
    End Sub

    Public Overloads Shared Sub RGBToGammaRGB(ByVal sRGB() As Single, ByVal sGammaRGB() As Single)
        RGBToGammaRGB(sRGB, sGammaRGB, 255)
    End Sub

    Public Overloads Shared Sub RGBToGammaRGB(ByVal sRGB() As Single, ByVal sGammaRGB() As Single, ByVal iScale As Integer)
        'Gamma correction. Input RGB in [0.0 1.0], output RGB [0 iScale].
        sGammaRGB(RGB.red) = GammaCorrection(sRGB(RGB.red)) * iScale
        sGammaRGB(RGB.green) = GammaCorrection(sRGB(RGB.green)) * iScale
        sGammaRGB(RGB.blue) = GammaCorrection(sRGB(RGB.blue)) * iScale
    End Sub

    Public Overloads Shared Sub GammaRGBToRGB(ByVal iGammaRGB() As Integer, ByVal sRGB() As Single, ByVal iScale As Integer)
        'Inverse gamma correction, single precision input.
        'Output RGB in [0.0 1.0], input RGB [0.0 iScale].
        'Input maybe same array as output.
        sRGB(RGB.red) = InverseGammaCorrection(CSng(iGammaRGB(RGB.red) / iScale))
        sRGB(RGB.green) = InverseGammaCorrection(CSng(iGammaRGB(RGB.green) / iScale))
        sRGB(RGB.blue) = InverseGammaCorrection(CSng(iGammaRGB(RGB.blue) / iScale))
    End Sub

    Public Overloads Shared Sub GammaRGBToRGB(ByVal iGammaRGB() As Integer, ByVal sRGB() As Single)
        GammaRGBToRGB(iGammaRGB, sRGB, 255)
    End Sub

    Public Overloads Shared Sub GammaRGBToRGB(ByVal sGammaRGB() As Single, ByVal sRGB() As Single, ByVal iScale As Integer)
        'Inverse gamma correction, single precision input.
        'Output RGB in [0.0 1.0], input RGB [0.0 iScale].
        'Input maybe same array as output.
        sRGB(RGB.red) = InverseGammaCorrection(sGammaRGB(RGB.red) / iScale)
        sRGB(RGB.green) = InverseGammaCorrection(sGammaRGB(RGB.green) / iScale)
        sRGB(RGB.blue) = InverseGammaCorrection(sGammaRGB(RGB.blue) / iScale)
    End Sub

    Public Overloads Shared Sub GammaRGBToRGB(ByVal sGammaRGB() As Single, ByVal sRGB() As Single)
        GammaRGBToRGB(sGammaRGB, sRGB, 255)
    End Sub

    Public Overloads Shared Sub IntegerRGBToRGB(ByVal iRGB() As Integer, ByVal sRGB() As Single)
        IntegerRGBToRGB(iRGB, sRGB, 255)
    End Sub

    Public Overloads Shared Sub IntegerRGBToRGB(ByVal iRGB() As Integer, ByVal sRGB() As Single, ByVal iScale As Integer)
        'Output RGB in [0.0 1.0], input RGB [0 iScale]. Simple scaling 
        sRGB(RGB.red) = CSng(iRGB(RGB.red) / iScale)
        sRGB(RGB.green) = CSng(iRGB(RGB.green) / iScale)
        sRGB(RGB.blue) = CSng(iRGB(RGB.blue) / iScale)
    End Sub

    Public Overloads Shared Sub RGBToIntegerRGB(ByVal sRGB() As Single, ByVal iRGB() As Integer)
        RGBToIntegerRGB(sRGB, iRGB, 255)
    End Sub

    Public Overloads Shared Sub RGBToIntegerRGB(ByVal sRGB() As Single, ByVal iRGB() As Integer, ByVal iScale As Integer)
        'Input RGB in [0.0 1.0], output RGB [0 iScale]. Simple scaling 
        iRGB(RGB.red) = CInt(Math.Round(sRGB(RGB.red) * iScale))
        iRGB(RGB.green) = CInt(Math.Round(sRGB(RGB.green) * iScale))
        iRGB(RGB.blue) = CInt(Math.Round(sRGB(RGB.blue) * iScale))
    End Sub

    Public Overloads Shared Sub RGBToLUTRGB(ByVal iRGBIn() As Integer, ByVal iRGBOut() As Integer, ByVal LUT(,) As Integer)
        'Apply LUT to integer from RGB or R*G*B* to the SAME color space. Input and output array maybe the same
        iRGBOut(RGB.red) = LUT(iRGBIn(RGB.red), RGB.red)
        iRGBOut(RGB.green) = LUT(iRGBIn(RGB.green), RGB.green)
        iRGBOut(RGB.blue) = LUT(iRGBIn(RGB.blue), RGB.blue)
    End Sub

    Public Overloads Shared Sub RGBToLUTRGB(ByVal iRGBIn() As Integer, ByVal sRGBOut() As Single, ByVal LUT(,) As Single)
        'Apply LUT to integer from RGB or R*G*B* to RGB. LUT must incorporate inverse gamma! Input and output array maybe the same
        sRGBOut(RGB.red) = LUT(iRGBIn(RGB.red), RGB.red)
        sRGBOut(RGB.green) = LUT(iRGBIn(RGB.green), RGB.green)
        sRGBOut(RGB.blue) = LUT(iRGBIn(RGB.blue), RGB.blue)
    End Sub

    Public Overloads Shared Function ToString(ByVal sColor() As Single) As String
        Return ToString(sColor, 3)
    End Function

    Public Overloads Shared Function ToString(ByVal sColor() As Single, _
                                       ByVal iPrecision As Integer) As String
        'Transform a 3 element color to a string for screen output
        Dim strFormat As String, i As Integer

        If iPrecision > 0 Then
            strFormat = "##0."
            For i = 1 To iPrecision
                strFormat = strFormat & "#"
            Next i
        Else
            strFormat = "#"
        End If
        Return sColor(0).ToString(strFormat) & " " & sColor(1).ToString(strFormat) & " " & sColor(2).ToString(strFormat)
    End Function

    Public Overloads Shared Function ToString(ByVal byColor() As Integer) As String
        'Transform a 3 element color to a string for screen output
        Return byColor(0) & " " & byColor(1) & " " & byColor(2)
    End Function

    Public Overloads Shared Function ToString(ByVal byColor() As Byte) As String
        'Transform a 3 element color to a string for screen output
        Return byColor(0) & " " & byColor(1) & " " & byColor(2)
    End Function

    Public Shared Sub ComputeColorSpaceTransform(ByVal sC1(,) As Single, _
                                                 ByVal sC2(,) As Single, _
                                                 ByRef sTransform(,) As Single, _
                                                 ByVal iTransformType As ColorSpaceTransform)
        'Compute the polynomial transforms from R3 to R3, both representing LINEAR color spaces.
        'This is based on a set of N color triplets in sC1, which have to be mapped to N color triplets in sC2.
        'This means both sC1 and sC2 are N x 3 matrices. 
        'This also means there are actually 3 transforms to be computed, i.e. one for every color coordinate.
        'The passed on iTransformType is equal to the number of terms in the transform (3 for linear,etc..), 
        'and can be termed the order of the transform.
        'The resulting polynomial transform  is stored in a 3 x M matrix, which can be right multiplied with
        'a color N-tuple to obtain the desired color tristimulus value.
        'NOTE: Usually the color transform is used to transform from RGB to sRGB or XYZ
        Dim iNrColors As Integer

        iNrColors = sC1.GetUpperBound(0) + 1
        Debug.WriteLine("Computing " & iTransformType & _
          "-term transform using " & iNrColors & " tristimulus values")

        'Compute sA using the input color triplets
        Dim sA(,) As Single = Nothing
        ColorTripletToNTuple(sC1, sA, iTransformType)
        Debug.WriteLine("Design matrix is " & vbNewLine & Algebra.ToString(sA))

        'Solve the 3 sets of linear equations (they are linear in their unknowns, even if they use
        'coefficient which are non-linear functions of the input color triplets!)
        'The resulting matrix needs to be transposed to be consistent with the fixed linear transform
        'matrices already used in the color class (they all use Tf.X = Y, NOT X.Tf = Y)
        Dim sX(iTransformType - 1, 2) As Single
        Algebra.Solve(sA, sX, sC2)

        ReDim sTransform(2, iTransformType - 1)
        Algebra.Transpose(sX, sTransform)
    End Sub

    Public Sub ImproveRGBTosRGBTransform(ByVal sInputColors(,) As Single, _
                                         ByVal sAimColors(,) As Single, _
                                         ByRef sTransform(,) As Single, _
                                         ByVal iMetricPower As Integer, _
                                         ByVal iNrRuns As Integer)
        'Improve the RGB to sRGB color space transform using non-linear optimization
        'IF necessary create other function to improve other types of transforms ....
        'Because we use delegates we need to set some global variables ...

        'The transform has 3 x iNrTerms parameters (or dimensions)
        'as the R,G and B transforms are no longer computed
        'independently as in the SVD method. The coefficients
        'must be in a big vector scoef which is filled before each
        'run
        Dim iNrColors As Integer = sInputColors.GetUpperBound(0) + 1
        Dim iNrTerms As Integer = sTransform.GetUpperBound(1) + 1
        Dim iNrDims As Integer = 3 * iNrTerms
        Dim i As Integer

        Debug.WriteLine("Improving " & iNrTerms & _
          "-term transform using " & iNrColors & " tristimulus values")

        iColorSpaceTransformMetricPower = iMetricPower
        Debug.WriteLine("Using " & iColorSpaceTransformMetricPower & "-metric in CIE Lab space")

        'Copy the colors. Transform them to CIELab to speed up things
        ReDim sColorSpaceTransformAimCIELabColors(sAimColors.GetUpperBound(0), sAimColors.GetUpperBound(1))
        For i = 0 To iNrColors - 1
            Dim ssRGB(2), sLab(2) As Single
            Algebra.GetMatrixRow(sAimColors, ssRGB, i)
            sRGBToCIELab(ssRGB, sLab)
            Algebra.SetMatrixRow(sColorSpaceTransformAimCIELabColors, sLab, i)
        Next

        'Compute the designmatrix using the input color triplets
        ColorTripletToNTuple(sInputColors, sColorSpaceTransformDesignMatrix, iNrTerms)
        'debug.writeline("Design matrix is " & vbNewLine & Algebra.ToString(sColorSpaceTransformDesignMatrix))

        'Initialize a simplex using the initial coefficient values
        'This simplex kind of defines the initial search space, so
        'this is important. One of the vertices is the SVD solution.
        'One coefficient of the other vertices is perturbed by
        'a multiplication by a factor between 0.5 and 1.5
        Dim iNrVertices As Integer = iNrDims + 1
        Dim sSimplex(iNrVertices - 1, iNrDims - 1) As Single
        Dim sSimplexVal(iNrVertices - 1) As Single
        Dim sCoef(iNrDims - 1) As Single
        Dim sBestCoef(iNrDims - 1) As Single
        Dim sBestSimplexVal As Single = 10000000
        Dim iRun, iBestRun As Integer
        Randomize()

        For iRun = 0 To iNrRuns - 1
            For i = 0 To iNrVertices - 1
                Algebra.MatrixToVector(sTransform, sCoef)
                If i > 0 Then
                    If Rnd() > 0.5 Then
                        sCoef(i - 1) = CSng(sCoef(i - 1) * (Rnd() * 0.75 + 0.25))
                    Else
                        sCoef(i - 1) = sCoef(i - 1) * (Rnd() * 4 + 1)
                    End If
                End If

                'Put in simplex and compute function value
                Algebra.SetMatrixRow(sSimplex, sCoef, i)
                sSimplexVal(i) = RGBTosRGBTransformError(sCoef)
            Next i

            If iRun = 0 Then
                Debug.WriteLine("Initial transform: " & Algebra.ToString(sTransform))
                Debug.WriteLine("Error: " & sSimplexVal(0))
            End If

            'Optimize
            Debug.WriteLine("Run " & iRun)
            Debug.WriteLine("Simplex errors: " & Algebra.ToString(sSimplexVal))
            Algebra.SolveNonLinear(sSimplex, sSimplexVal, 32768, AddressOf Me.RGBTosRGBTransformError)

            If sSimplexVal(0) < sBestSimplexVal Then
                sBestSimplexVal = sSimplexVal(0)
                Algebra.GetMatrixRow(sSimplex, sBestCoef, 0)
                iBestRun = iRun
            End If
        Next

        'Copy result back to transform
        Debug.WriteLine("Best run is " & iBestRun & " , with error " & sBestSimplexVal)
        Algebra.VectorToMatrix(sBestCoef, sTransform)
    End Sub

    Public Function RGBTosRGBTransformError(ByVal sCoef() As Single) As Single
        'Return the error to be minimized. Uses global variables!!!!!
        'Any value for the errormetric above 3 means max-metric!
        Dim iNrColors As Integer = ColorUtil.sColorSpaceTransformAimCIELabColors.GetUpperBound(0) + 1
        Dim iNrTerms As Integer = ColorUtil.sColorSpaceTransformDesignMatrix.GetUpperBound(1) + 1
        Dim sError As Single

        Dim i As Integer, j As Integer, k As Integer
        Dim iCoefNr As Integer
        Dim ssRGB(2), sLab(2), sAimLab(2) As Single

        sError = 0
        For i = 0 To iNrColors - 1
            iCoefNr = 0
            For j = 0 To 2
                sAimLab(j) = ColorUtil.sColorSpaceTransformAimCIELabColors(i, j)
                ssRGB(j) = 0
                For k = 0 To iNrTerms - 1
                    ssRGB(j) = ssRGB(j) + sCoef(iCoefNr) * ColorUtil.sColorSpaceTransformDesignMatrix(i, k)
                    iCoefNr = iCoefNr + 1
                Next k
            Next j

            sRGBToCIELab(ssRGB, sLab)
            If iColorSpaceTransformMetricPower < 4 Then
                sError = CSng(sError + Math.Pow(CIELabTodE(sAimLab, sLab), iColorSpaceTransformMetricPower))
            Else
                sError = Math.Max(sError, CIELabTodE(sAimLab, sLab))
            End If
        Next i

        If iColorSpaceTransformMetricPower < 4 Then
            sError = CSng(Math.Pow(sError / iNrColors, 1 / iColorSpaceTransformMetricPower))
        End If
        Return sError
    End Function

    Private Overloads Shared Sub ColorTripletToNTuple(ByVal sCList(,) As Single, _
                                                      ByRef sCNLList(,) As Single, _
                                                      ByVal iNrTerms As Integer)
        'Computes N-tuples for several color triplets at once.
        Dim iNrColors As Integer = sCList.GetUpperBound(0) + 1
        Dim sC(2), sCNL(iNrTerms - 1) As Single, i As Integer
        ReDim sCNLList(iNrColors - 1, iNrTerms - 1)

        For i = 0 To iNrColors - 1
            Algebra.GetMatrixRow(sCList, sC, i)
            ColorTripletToNTuple(sC, sCNL, iNrTerms)

            'Copy color N-tuple to sA
            Algebra.SetMatrixRow(sCNLList, sCNL, i)
        Next
    End Sub

    Private Overloads Shared Sub ColorTripletToNTuple(ByVal sC() As Single, _
                                                      ByVal sCNL() As Single, _
                                                      ByVal iNrTerms As Integer)
        'Compute a color N-tuple using higher order and cross terms of the input color triplet.
        'N depends on the number of elements in sCNL.
        'This N-tuple can then be multiplied by a transform matrix representing the corresponding 
        'polynomial transform to obtain a new color triplet.

        'All base function sets are at least linear in RGB
        Array.Copy(sC, sCNL, 3)

        Select Case iNrTerms
            Case ColorSpaceTransform.Linear
                'Linear transform
            Case ColorSpaceTransform.NonLinear20Term
                ' 20 terms:
                ' tf = a0 + a1*X + a2*y + a3*z + a4*X*y + a5*y*z + a6*z*X + a7*X^2 + a8*y^2 + a9*z^2 + a10 * X*y*z
                '      + a11*X^3 + a12*y^3 + a13*z^3 + a14*X*y^2 + a15*X^2*y + a16*y*z^2 + a17*y^2*z + a18*z*X^2
                '      + a19*z^2*X
                sCNL(3) = sCNL(0) * sCNL(1)
                sCNL(4) = sCNL(2) * sCNL(1)
                sCNL(5) = sCNL(0) * sCNL(2)

                sCNL(6) = sCNL(0) * sCNL(0)
                sCNL(7) = sCNL(1) * sCNL(1)
                sCNL(8) = sCNL(2) * sCNL(2)

                sCNL(9) = sCNL(0) * sCNL(1) * sCNL(2)
                sCNL(10) = sCNL(6) * sCNL(0)
                sCNL(11) = sCNL(7) * sCNL(1)
                sCNL(12) = sCNL(8) * sCNL(2)

                sCNL(13) = sCNL(3) * sCNL(1)
                sCNL(14) = sCNL(3) * sCNL(0)
                sCNL(15) = sCNL(4) * sCNL(2)
                sCNL(16) = sCNL(4) * sCNL(1)
                sCNL(17) = sCNL(5) * sCNL(0)
                sCNL(18) = sCNL(5) * sCNL(2)
                sCNL(19) = 1.0#
            Case ColorSpaceTransform.NonLinear14Term
                ' 14 terms:
                ' tf = a13 + a0*X + a1*y + a2*z + a3*X*y + a4*y*z + a5*z*X +
                '      a6*X^2 + a7*y^2 + a8*z^2 + a9 * X*y*z
                '      + a10*X^3 + a11*y^3 + a12*z^3
                sCNL(3) = sCNL(0) * sCNL(1)
                sCNL(4) = sCNL(2) * sCNL(1)
                sCNL(5) = sCNL(0) * sCNL(2)

                sCNL(6) = sCNL(0) * sCNL(0)
                sCNL(7) = sCNL(1) * sCNL(1)
                sCNL(8) = sCNL(2) * sCNL(2)

                sCNL(9) = sCNL(0) * sCNL(1) * sCNL(2)
                sCNL(10) = sCNL(6) * sCNL(0)
                sCNL(11) = sCNL(7) * sCNL(1)
                sCNL(12) = sCNL(8) * sCNL(2)

                sCNL(13) = 1.0
            Case ColorSpaceTransform.NonLinear11Term
                ' 11 terms:
                ' tf = a10 + a0*X + a1*y + a2*z + a3*X*y + a4*y*z + a5*z*X +
                '      a6*X^2 + a7*y^2 + a8*z^2 + a9*X*y*z
                sCNL(3) = sCNL(0) * sCNL(1)
                sCNL(4) = sCNL(2) * sCNL(1)
                sCNL(5) = sCNL(0) * sCNL(2)

                sCNL(6) = sCNL(0) * sCNL(0)
                sCNL(7) = sCNL(1) * sCNL(1)
                sCNL(8) = sCNL(2) * sCNL(2)

                sCNL(9) = sCNL(0) * sCNL(1) * sCNL(2)
                sCNL(10) = 1.0#
            Case ColorSpaceTransform.NonLinear6Term
                ' 6 terms:
                ' tf = a0*X + a1*y + a2*z + a3*X*y + a4*y*z + a5*z*X
                sCNL(3) = sCNL(0) * sCNL(1)
                sCNL(4) = sCNL(2) * sCNL(1)
                sCNL(5) = sCNL(0) * sCNL(2)
            Case ColorSpaceTransform.NonLinear8Term
                ' 8 terms:
                ' tf = a7 + a0*X + a1*y + a2*z + a3*X*y +
                '      a4*y*z + a5*z*X + a6*X*y*z
                sCNL(3) = sCNL(0) * sCNL(1)
                sCNL(4) = sCNL(2) * sCNL(1)
                sCNL(5) = sCNL(0) * sCNL(2)

                sCNL(6) = sCNL(0) * sCNL(1) * sCNL(2)
                sCNL(7) = 1.0
            Case ColorSpaceTransform.NonLinear9Term
                ' 9 terms:
                ' tf = a0*X + a1*y + a2*z + a3*X*y + a4*y*z + a5*z*X +
                '      a6*X^2 + a7*y^2 + a8*z^2
                sCNL(3) = sCNL(0) * sCNL(1)
                sCNL(4) = sCNL(2) * sCNL(1)
                sCNL(5) = sCNL(0) * sCNL(2)

                sCNL(6) = sCNL(0) * sCNL(0)
                sCNL(7) = sCNL(1) * sCNL(1)
                sCNL(8) = sCNL(2) * sCNL(2)
            Case Else
        End Select

        'debug.writeline("Converted color triplet:" & _
        '  vbNewLine & Algebra.ToString(sC) & _
        '  "to " & _
        '  vbNewLine & Algebra.ToString(sCNL))
    End Sub

    '------------------- Some functions for use with Windows GDI+ Get and SetPixel ----------------------

    Public Overloads Shared Function GammaRGBToColor(ByVal sGammaRGB() As Single) As Color
        'Return a color structure from gamma-corrected RGB
        Dim objColor As Color
        objColor = Color.FromArgb(CInt(sGammaRGB(RGB.red)), CInt(sGammaRGB(RGB.green)), CInt(sGammaRGB(RGB.blue)))
    End Function

    Public Overloads Shared Function GammaRGBToColor(ByVal byGammaRGB() As Byte) As Color
        'Return a color structure from gamma-corrected RGB
        Return Color.FromArgb(byGammaRGB(RGB.red), byGammaRGB(RGB.green), byGammaRGB(RGB.blue))
    End Function

    Public Shared Function RGBToColor(ByVal sRGB() As Single) As Color
        'Return a color structure from LINEAR RGB
        Dim iRGB(2) As Integer
        RGBToGammaRGB(sRGB, iRGB)
        Return Color.FromArgb(iRGB(RGB.red), iRGB(RGB.green), iRGB(RGB.blue))
    End Function

    Public Overloads Shared Sub ColorToRGB(ByVal objColor As Color, ByVal sRGB() As Single)
        sRGB(RGB.red) = objColor.R
        sRGB(RGB.green) = objColor.G
        sRGB(RGB.blue) = objColor.B
        GammaRGBToRGB(sRGB, sRGB)
    End Sub

    Public Overloads Shared Function RGBColorTosRGBColor(ByVal objColor As Color, ByVal sTransform(,) As Single) As Color
        'Conveniece function that tranforms an color struct containing gamma-corrected
        'RGB values to sRGB values
        Dim sRGB(2), sSRGB(2) As Single

        ColorToRGB(objColor, sRGB)
        TransformTristimulus(sRGB, sSRGB, sTransform)
        Return RGBToColor(sSRGB)
    End Function
End Class
